Vittorio Amoruso
Nicola Zucchia
Erion Islamay
R, Rmd, HTML, CSS
library(boot)
library(car)
library(ellipse)
library(ggplot2)
library(gridExtra)
library(corrplot)
library(RColorBrewer)
library(GGally)
library(cluster)
Text Edit
data presente all’interno del file in formato testuale
data conforme alla convenzione csv
Header
mydata<-read.csv("imports-85.data",header = FALSE)
Correttezza e Struttura
dim(mydata)
## [1] 205 26
Il Dataframe presenta 205 osservazioni in 26 variabili
Visualizzazione del Contenuto
# View(ds)
head(mydata, 5)
## V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13
## 1 3 ? alfa-romero gas std two convertible rwd front 88.6 168.8 64.1 48.8
## 2 3 ? alfa-romero gas std two convertible rwd front 88.6 168.8 64.1 48.8
## 3 1 ? alfa-romero gas std two hatchback rwd front 94.5 171.2 65.5 52.4
## 4 2 164 audi gas std four sedan fwd front 99.8 176.6 66.2 54.3
## 5 2 164 audi gas std four sedan 4wd front 99.4 176.6 66.4 54.3
## V14 V15 V16 V17 V18 V19 V20 V21 V22 V23 V24 V25 V26
## 1 2548 dohc four 130 mpfi 3.47 2.68 9 111 5000 21 27 13495
## 2 2548 dohc four 130 mpfi 3.47 2.68 9 111 5000 21 27 16500
## 3 2823 ohcv six 152 mpfi 2.68 3.47 9 154 5000 19 26 16500
## 4 2337 ohc four 109 mpfi 3.19 3.40 10 102 5500 24 30 13950
## 5 2824 ohc five 136 mpfi 3.19 3.40 8 115 5500 18 22 17450
Assegnazione dei nomi - CamelCase
names(mydata) <- c('Symboling', 'NormalizedLosses', 'Make', 'FuelType', 'Aspiration', 'NumOfDoors', 'BodyStyle', 'DriveWheels', 'EngineLocation', 'WheelBase', 'Length', 'Width', 'Height', 'CurbWeight', 'EngineType', 'NumOfCylinders', 'EngineSize', 'FuelSystem', 'Bore', 'Stroke', 'CompressionRatio', 'Horsepower', 'PeakRpm', 'CityMpg', 'HighwayMpg', 'Price')
Qualitative vs Quantitative
str(mydata)
## 'data.frame': 205 obs. of 26 variables:
## $ Symboling : int 3 3 1 2 2 2 1 1 1 0 ...
## $ NormalizedLosses: chr "?" "?" "?" "164" ...
## $ Make : chr "alfa-romero" "alfa-romero" "alfa-romero" "audi" ...
## $ FuelType : chr "gas" "gas" "gas" "gas" ...
## $ Aspiration : chr "std" "std" "std" "std" ...
## $ NumOfDoors : chr "two" "two" "two" "four" ...
## $ BodyStyle : chr "convertible" "convertible" "hatchback" "sedan" ...
## $ DriveWheels : chr "rwd" "rwd" "rwd" "fwd" ...
## $ EngineLocation : chr "front" "front" "front" "front" ...
## $ WheelBase : num 88.6 88.6 94.5 99.8 99.4 ...
## $ Length : num 169 169 171 177 177 ...
## $ Width : num 64.1 64.1 65.5 66.2 66.4 66.3 71.4 71.4 71.4 67.9 ...
## $ Height : num 48.8 48.8 52.4 54.3 54.3 53.1 55.7 55.7 55.9 52 ...
## $ CurbWeight : int 2548 2548 2823 2337 2824 2507 2844 2954 3086 3053 ...
## $ EngineType : chr "dohc" "dohc" "ohcv" "ohc" ...
## $ NumOfCylinders : chr "four" "four" "six" "four" ...
## $ EngineSize : int 130 130 152 109 136 136 136 136 131 131 ...
## $ FuelSystem : chr "mpfi" "mpfi" "mpfi" "mpfi" ...
## $ Bore : chr "3.47" "3.47" "2.68" "3.19" ...
## $ Stroke : chr "2.68" "2.68" "3.47" "3.40" ...
## $ CompressionRatio: num 9 9 9 10 8 8.5 8.5 8.5 8.3 7 ...
## $ Horsepower : chr "111" "111" "154" "102" ...
## $ PeakRpm : chr "5000" "5000" "5000" "5500" ...
## $ CityMpg : int 21 21 19 24 18 19 19 19 17 16 ...
## $ HighwayMpg : int 27 27 26 30 22 25 25 25 20 22 ...
## $ Price : chr "13495" "16500" "16500" "13950" ...
mydata$Symboling <- factor(mydata$Symboling)
mydata$Make <- factor(mydata$Make)
mydata$FuelType <- factor(mydata$FuelType)
mydata$Aspiration <- factor(mydata$Aspiration)
mydata$NumOfDoors <- factor(mydata$NumOfDoors)
mydata$BodyStyle <- factor(mydata$BodyStyle)
mydata$DriveWheels <- factor(mydata$DriveWheels)
mydata$EngineLocation <- factor(mydata$EngineLocation)
mydata$EngineType <- factor(mydata$EngineType)
mydata$NumOfCylinders <- factor(mydata$NumOfCylinders)
mydata$FuelSystem <- factor(mydata$FuelSystem)
mydata$NormalizedLosses <- as.numeric(mydata$NormalizedLosses)
mydata$Bore <- as.numeric(mydata$Bore)
mydata$Stroke <- as.numeric(mydata$Stroke)
mydata$Horsepower <- as.numeric(mydata$Horsepower)
mydata$PeakRpm <- as.numeric(mydata$PeakRpm)
mydata$Price <- as.numeric(mydata$Price)
str(mydata)
## 'data.frame': 205 obs. of 26 variables:
## $ Symboling : Factor w/ 6 levels "-2","-1","0",..: 6 6 4 5 5 5 4 4 4 3 ...
## $ NormalizedLosses: num NA NA NA 164 164 NA 158 NA 158 NA ...
## $ Make : Factor w/ 22 levels "alfa-romero",..: 1 1 1 2 2 2 2 2 2 2 ...
## $ FuelType : Factor w/ 2 levels "diesel","gas": 2 2 2 2 2 2 2 2 2 2 ...
## $ Aspiration : Factor w/ 2 levels "std","turbo": 1 1 1 1 1 1 1 1 2 2 ...
## $ NumOfDoors : Factor w/ 3 levels "?","four","two": 3 3 3 2 2 3 2 2 2 3 ...
## $ BodyStyle : Factor w/ 5 levels "convertible",..: 1 1 3 4 4 4 4 5 4 3 ...
## $ DriveWheels : Factor w/ 3 levels "4wd","fwd","rwd": 3 3 3 2 1 2 2 2 2 1 ...
## $ EngineLocation : Factor w/ 2 levels "front","rear": 1 1 1 1 1 1 1 1 1 1 ...
## $ WheelBase : num 88.6 88.6 94.5 99.8 99.4 ...
## $ Length : num 169 169 171 177 177 ...
## $ Width : num 64.1 64.1 65.5 66.2 66.4 66.3 71.4 71.4 71.4 67.9 ...
## $ Height : num 48.8 48.8 52.4 54.3 54.3 53.1 55.7 55.7 55.9 52 ...
## $ CurbWeight : int 2548 2548 2823 2337 2824 2507 2844 2954 3086 3053 ...
## $ EngineType : Factor w/ 7 levels "dohc","dohcv",..: 1 1 6 4 4 4 4 4 4 4 ...
## $ NumOfCylinders : Factor w/ 7 levels "eight","five",..: 3 3 4 3 2 2 2 2 2 2 ...
## $ EngineSize : int 130 130 152 109 136 136 136 136 131 131 ...
## $ FuelSystem : Factor w/ 8 levels "1bbl","2bbl",..: 6 6 6 6 6 6 6 6 6 6 ...
## $ Bore : num 3.47 3.47 2.68 3.19 3.19 3.19 3.19 3.19 3.13 3.13 ...
## $ Stroke : num 2.68 2.68 3.47 3.4 3.4 3.4 3.4 3.4 3.4 3.4 ...
## $ CompressionRatio: num 9 9 9 10 8 8.5 8.5 8.5 8.3 7 ...
## $ Horsepower : num 111 111 154 102 115 110 110 110 140 160 ...
## $ PeakRpm : num 5000 5000 5000 5500 5500 5500 5500 5500 5500 5500 ...
## $ CityMpg : int 21 21 19 24 18 19 19 19 17 16 ...
## $ HighwayMpg : int 27 27 26 30 22 25 25 25 20 22 ...
## $ Price : num 13495 16500 16500 13950 17450 ...
Eliminazione V. non di Interesse
mydata <- mydata[ , - c(1, 2) ]
Gestione degli NA, ‘?’
mydata <- na.omit(mydata)
cleaner <- function(ds, sc, show = FALSE) {
# ds = DataSet , sc = SpecialCharacter
unlist <- as.numeric( )
for (i in 1:dim(ds)[1])
{
dr <- ds[ i , ]
if ( any(dr == sc) ) unlist[length(unlist) + 1] <- i
}
if (show == TRUE)
{
print('Le unita eliminate sono.. ')
print(unlist)
}
return(ds[ - unlist , ])
}
mydata <- cleaner(mydata, '?', show = TRUE)
## [1] "Le unita eliminate sono.. "
## [1] 27 57
Make <- table(mydata$Make)
FuelType <- table(mydata$FuelType)
Aspiration <- table(mydata$Aspiration)
NumOfDoors <- table(mydata$NumOfDoors)
BodyStyle <- table(mydata$BodyStyle)
DriveWheels <- table(mydata$DriveWheels)
EngineLocation <- table(mydata$EngineLocation)
EngineType <- table(mydata$EngineType)
NumOfCylinders <- table(mydata$NumOfCylinders)
FuelSystem <- table(mydata$FuelSystem)
barplot(BodyStyle, col = c("grey", "grey", "lightblue", "#eb8060", "grey"), border=NA, main = "BodyStyle")
barplot(DriveWheels, col = c("grey", "#eb8060", "lightblue"), border=NA, main = "DriveWheels")
barplot(EngineType, col = c("grey", "grey", "grey", "#eb8060", "grey", "grey", "grey"), border=NA, main = "EngineType")
barplot(NumOfCylinders, col = c("grey", "grey", "#eb8060", "grey", "grey", "grey", "grey"), border=NA, main = "NumOfCylinders")
barplot(FuelSystem, col = c("grey", "lightblue", "grey", "grey", "grey", "#eb8060", "grey", "grey"), border=NA, main = "FuelSystem")
table(mydata$FuelType, mydata$Aspiration)
##
## std turbo
## diesel 6 13
## gas 152 22
chisq.test(table(mydata$FuelType, mydata$Aspiration))
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table(mydata$FuelType, mydata$Aspiration)
## X-squared = 32.238, df = 1, p-value = 1.364e-08
table(mydata$DriveWheels, mydata$FuelType)
##
## diesel gas
## 4wd 0 8
## fwd 8 106
## rwd 11 60
chisq.test(table(mydata$DriveWheels, mydata$FuelType))
##
## Pearson's Chi-squared test
##
## data: table(mydata$DriveWheels, mydata$FuelType)
## X-squared = 4.4523, df = 2, p-value = 0.1079
#plot maximum price according to the automobile maker i.e., brand
maxPrice <- aggregate(mydata$Price , by = list(mydata$Make), FUN = "max")
names(maxPrice) <- c("Brand", "Price")
maxPrice <- maxPrice[ order(maxPrice$Price) , ]
maxPrice$Brand <- factor(maxPrice$Brand, levels = maxPrice$Brand)
ggplot(data = maxPrice, aes(x = Price, y = Brand)) +
geom_bar(stat = "identity", fill = "lightblue") +
geom_text(aes(label = paste0("$ ", Price)), hjust = -0.05) +
coord_cartesian(xlim = c(0, 52000)) +
labs(title = "Max price of automobile in each Brand")
meanPrice <- aggregate(mydata$Price, by = list(mydata$Make), FUN = "mean")
names(meanPrice) <- c("Brand","Price")
meanPrice <- meanPrice[ order(meanPrice$Price) , ]
meanPrice$Brand <- factor(meanPrice$Brand, levels = meanPrice$Brand)
ggplot(data = meanPrice, aes(x = Price, y = Brand)) +
geom_bar(stat = "identity", fill = "Lightblue") +
geom_text(aes(label = paste0("$ ", round(Price))), hjust = -0.05) +
coord_cartesian(xlim = c(0,48000)) +
labs(title = "Mean price of automobile in each Brand")
mydata.num<- mydata[,c("WheelBase","Length","Width","Height","CurbWeight",
"EngineSize","Bore","Stroke","CompressionRatio","Horsepower",
"PeakRpm","CityMpg","HighwayMpg","Price")]
summary(mydata.num)
## WheelBase Length Width Height
## Min. : 86.60 Min. :141.1 Min. :60.30 Min. :47.80
## 1st Qu.: 94.50 1st Qu.:166.3 1st Qu.:64.10 1st Qu.:52.00
## Median : 97.00 Median :173.2 Median :65.40 Median :54.10
## Mean : 98.92 Mean :174.3 Mean :65.89 Mean :53.87
## 3rd Qu.:102.40 3rd Qu.:184.6 3rd Qu.:66.90 3rd Qu.:55.70
## Max. :120.90 Max. :208.1 Max. :72.00 Max. :59.80
## CurbWeight EngineSize Bore Stroke
## Min. :1488 Min. : 61.0 Min. :2.540 Min. :2.070
## 1st Qu.:2145 1st Qu.: 98.0 1st Qu.:3.150 1st Qu.:3.110
## Median :2414 Median :120.0 Median :3.310 Median :3.290
## Mean :2562 Mean :128.1 Mean :3.331 Mean :3.249
## 3rd Qu.:2952 3rd Qu.:146.0 3rd Qu.:3.590 3rd Qu.:3.410
## Max. :4066 Max. :326.0 Max. :3.940 Max. :4.170
## CompressionRatio Horsepower PeakRpm CityMpg
## Min. : 7.00 Min. : 48.0 Min. :4150 Min. :13.00
## 1st Qu.: 8.50 1st Qu.: 70.0 1st Qu.:4800 1st Qu.:19.00
## Median : 9.00 Median : 95.0 Median :5100 Median :25.00
## Mean :10.14 Mean :103.5 Mean :5100 Mean :25.33
## 3rd Qu.: 9.40 3rd Qu.:116.0 3rd Qu.:5500 3rd Qu.:30.00
## Max. :23.00 Max. :262.0 Max. :6600 Max. :49.00
## HighwayMpg Price
## Min. :16.00 Min. : 5118
## 1st Qu.:25.00 1st Qu.: 7738
## Median :30.00 Median :10245
## Mean :30.79 Mean :13285
## 3rd Qu.:34.00 3rd Qu.:16515
## Max. :54.00 Max. :45400
ggpairs(mydata.num,
columns = c("Length","Width","CurbWeight", "Price"),
aes(color = mydata$FuelType, # Color by group (cat. variable)
alpha = 0.5), upper = list(continuous = "points"))
ggpairs(mydata.num,
columns = c("EngineSize","Horsepower", "Price"),
aes(color = mydata$FuelType, # Color by group (cat. variable)
alpha = 0.5), upper = list(continuous = "points"))
ggpairs(mydata.num,
columns = c("CityMpg","HighwayMpg","Price"),
aes(color = mydata$FuelType, # Color by group (cat. variable)
alpha = 0.5), upper = list(continuous = "points"))
corr<- round(cor(mydata.num),2)
corrplot(corr, type="upper", order="hclust",
col=brewer.pal(n=8, name="RdYlBu"))
Studiamo le distribuzioni delle principali variabili quantitative interessanti mostrando qualche istogramma e curve di densità.
Si osserva che, per alcune variabili, la scarsa ampiezza campionaria non rende efficace la visualizzazione dei dati tramite istogrammi.
ggplot(mydata, aes(x = "", y = Price)) +
geom_boxplot(width = 0.4, colour="darkblue", fill="lightblue") +
geom_jitter( colour="red",
width = 0.1, size = 1) +
scale_color_manual(values = c("#00AFBB", "#E7B800")) +
labs(x = NULL)
ggplot(mydata, aes(x = "", y = log(Price))) +
geom_boxplot(width = 0.4, colour="darkblue", fill="lightblue") +
geom_jitter( colour="red",
width = 0.1, size = 1) +
scale_color_manual(values = c("#00AFBB", "#E7B800")) +
labs(x = NULL)
Come si nota dal boxplot, la variabile Price presenta una asimmetria marcata e molti outliers nella coda destra della distribuzione. Si osserva come fare una trasformazione logaritmica renda la distribuzione più regolare. Per esserne sicuri facciamo un confronto con la Gaussiana.
#si vede che il prezzo è fortemente assimetrico quindi forse da trasformare in log
#plot(density(mydata$price))
ggplot(mydata, aes(x=Price)) +
geom_histogram(aes(y=..density..), colour="darkblue", fill="lightblue")+
geom_density(alpha=.2, fill="#FF6666")
ggplot(mydata, aes(x=log(Price))) +
geom_histogram(aes(y=..density..), colour="darkblue", fill="lightblue")+
geom_density(alpha=.2, fill="#FF6666")
La variabile logPrice sembra avere una tendenza trimodale a partire dal grafico della sua densità.
Si osserva come fare una trasformazione logaritmica renda la distribuzione più regolare. Per esserne sicuri facciamo un confronto con la Gaussiana.
Testiamo la normalità della variabile price e della sua trasformata logPrice.
plot(ecdf(scale(mydata$Price)), cex=0.5, main="Price vs Pnorm")
curve(pnorm(x), add=TRUE)
plot(ecdf(scale(log(mydata$Price))), cex=0.5, main="logPrice vs Pnorm")
curve(pnorm(x), add=TRUE)
Il confronto con la funzione di ripartizione empirica suggerisce che una la variabile trasformata tenda a distribuirsi più normalmente della variabile di partenza.
qqnorm(mydata$Price, ylab = "Price")
qqline(mydata$Price)
qqnorm(log(mydata$Price), ylab = "LogPrice")
qqline(log(mydata$Price))
Facciamo un confronto anche di altre variabili quantitative che potrebbero invece avere una distribuzione normale
#lunghezza
qqnorm(mydata$Length, ylab = "Length")
qqline(mydata$Length)
#largezza
qqnorm(mydata$Width, ylab = "Width")
qqline(mydata$Width)
#altezza
qqnorm(mydata$Height, ylab = "Height")
qqline(mydata$Height)
A partire dai qqplot si deduce che le variabili Length e Height hanno una tendenza normale, mentre la variabile Width no.
Cominciamo a cercare un po di relazioni fra le variabili. Consideriamo il prezzo come VARIABILE RISPOSTA e cambiamo la VARIABILE COVARIATA di volta in volta.
confronto tra due insiemi di dati osservati:
ggplot(mydata, aes(x=Aspiration, y=log(Price), fill=Aspiration))+
geom_boxplot() + theme_gray(base_size = 14)
mean(log(mydata$Price[mydata$Aspiration=="std"]))
## [1] 9.287065
mean(log(mydata$Price[mydata$Aspiration=="turbo"]))
## [1] 9.642255
I boxplot mostrano che può esserci una correlazione significativa tra il tipo di aspirazione e il prezzo.
ggplot(mydata, aes(x=FuelType, y=log(Price), fill=FuelType))+
geom_boxplot()
#i diesel costao in media un po di piu
mean(log(mydata$Price[mydata$FuelType=="diesel"]))
## [1] 9.571661
mean(log(mydata$Price[mydata$FuelType=="gas"]))
## [1] 9.327435
Notiamo che ancora una volta c’è una differenza di medie tra le auto a benzina e quelle a diesel. Indaghiamo sulla significatività di tale differenza.
1) PRICE vs FUELTYPE
t.test(log(mydata$Price[mydata$FuelType=="diesel"]),
log(mydata$Price[mydata$FuelType=="gas"]), alternative = "two.sided")
##
## Welch Two Sample t-test
##
## data: log(mydata$Price[mydata$FuelType == "diesel"]) and log(mydata$Price[mydata$FuelType == "gas"])
## t = 2.0289, df = 22.314, p-value = 0.05457
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.005213846 0.493666882
## sample estimates:
## mean of x mean of y
## 9.571661 9.327435
Siccome il p-value è attorno a 0.05, siamo all’interno dell’intervallo di confidenza al 95%, e quindi non abbiamo prove sufficienti per rifiutare l’ipotesi nulla che il prezzo delle auto diesel non sia significativamente diverso da quello delle macchine a benzina.
Passiamo a vedere se la differenza per tipo di aspirazione è significativa
2) PRICE vs ASPIRATION
t.test(log(mydata$Price[mydata$Aspiration=="std"]),
log(mydata$Price[mydata$Aspiration=="turbo"]), alternative = "two.sided")
##
## Welch Two Sample t-test
##
## data: log(mydata$Price[mydata$Aspiration == "std"]) and log(mydata$Price[mydata$Aspiration == "turbo"])
## t = -4.6963, df = 65.777, p-value = 1.393e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.5062036 -0.2041761
## sample estimates:
## mean of x mean of y
## 9.287065 9.642255
Siccome il p-value è molto minore di 0.05, abbiamo prove sufficienti per rifiutare l’ipotesi nulla che il prezzo delle auto ad aspirazione standard non è significativamente diverso da quello delle macchine ad aspirazione turbo. In altre parole, il test ci dice che le macchina Turbo sono significativamente più costose delle macchine Standard.
3) PRICE vs BODYSTYLE
BodyStyle
##
## convertible hardtop hatchback sedan wagon
## 6 8 63 92 24
ggplot(mydata, aes(x=BodyStyle, y=log(Price), fill=BodyStyle))+
geom_boxplot()
#si vede che le macchine con hatcback (bagagliaio alto) sono piu economiche delle altre
Notando che le medie variano tra i vari stili di vetture, indaghiamo con il test ANOVA la significatività di tali differenze.
#ANOVA
df_aov <- aov(log(Price) ~ BodyStyle, data = mydata) #Fischer's classic ANOVA function
summary(df_aov)
## Df Sum Sq Mean Sq F value Pr(>F)
## BodyStyle 4 7.93 1.9820 8.822 1.51e-06 ***
## Residuals 188 42.24 0.2247
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(df_aov)
## Call:
## aov(formula = log(Price) ~ BodyStyle, data = mydata)
##
## Terms:
## BodyStyle Residuals
## Sum of Squares 7.92787 42.23654
## Deg. of Freedom 4 188
##
## Residual standard error: 0.4739857
## Estimated effects may be unbalanced
Dal test possiamo notare che c’è una differenza significativa, fra le medie dei gruppi. Dal boxplot possiamo assumere che almeno il prezzo delle auto hatcback sia sensibilmente differente dal prezzo delle altre auto.
4) PRICE vs DRIVE.WHEELS
DriveWheels
##
## 4wd fwd rwd
## 8 114 71
ggplot(mydata, aes(x=DriveWheels, y=log(Price), fill=DriveWheels))+
geom_boxplot()
#ANOVA
wheel_price_aov <- aov(log(Price) ~ DriveWheels, data = mydata)
summary(wheel_price_aov)
## Df Sum Sq Mean Sq F value Pr(>F)
## DriveWheels 2 24.12 12.060 87.98 <2e-16 ***
## Residuals 190 26.05 0.137
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(wheel_price_aov)
## Call:
## aov(formula = log(Price) ~ DriveWheels, data = mydata)
##
## Terms:
## DriveWheels Residuals
## Sum of Squares 24.11978 26.04462
## Deg. of Freedom 2 190
##
## Residual standard error: 0.3702391
## Estimated effects may be unbalanced
Il F-value è molto piccolo e questo porta a rifiutare l’ipotesi nulla che non vi sia sostanziale variabilità nel prezzo a seconda del tipo di trazione.
Vogliamo visualizzare alcune delle relazioni che intercorrono fra più di due variabili del dataset.
Vogliamo vedere la relazione che intercorre tra il prezzo delle auto e il loro consumo di carburante in città insieme alla tipologia di carburante.
Possiamo ipotizzare che le auto costose facciano meno chilometri con un litro, e che le machine diesel costino in media di più.
ggplot(mapping = aes(x=CityMpg, y=log(Price)), data=mydata) +
geom_point( alpha=0.4) + geom_smooth(aes(colour=FuelType),method = "lm", se=F) +
theme_gray(base_size = 14)
Come ipotizzavamo, esiste una correlazione negativa tra prezzo e consumo in città.
Adesso siamo interessati a vedere il collegamento tra prezzo, grandezza del motore e numero di cilindri.
Anche qui ci aspettiamo che ci sia una correlazione positiva: i motori grossi costano di più e hanno più cilindri.
ggplot(mapping = aes(x=EngineSize, y=log(Price)), data=mydata) +
geom_point( alpha=0.4) + geom_smooth(aes(colour=NumOfCylinders),method = "lm", se=F) +
theme_gray(base_size = 14)
Anche qui l’ipotesi sembra essere corretta.
Un’altra analisi interessante è vedere che relazione c’è tra prezzo e numero di cavalli. E anche qui cerchiamo di scoprire se cambia qualcosa a seconda del carburante.
# c'è una relazione interessante tra horsepower e prezzo come ci potevamo aspettare
ggplot(mapping = aes(x=Horsepower, y=log(Price)), data=mydata) +
geom_point(alpha=0.4) + geom_smooth(aes(colour=FuelType), method = "lm", se=F)
In sintesi notiamo che più cavalli abbiamo , più costa l’auto, ma quelle diesel all’aumentare dei cavalli sembrano essere più costose.
Insomma se voglio una macchina con tanti cavalli meglio prenderla a benzina.
Infine analizziamo come varia il prezzo a seconda del numero di cavalli e il tipo di aspirazione.
Anche qui ci aspettiamo che il prezzo sia più alto al crescere del numero dei cavalli e che l’aspirazione turbo abbia una tendenza ad essere più costosa.
ggplot(mapping = aes(x=Horsepower, y=log(Price)), data=mydata) +
geom_point(alpha=0.4) + geom_smooth(aes(colour=Aspiration), method = "lm", se=F)
E invece notiamo una cosa interessante: il prezzo dell’aspirazione turbo è più alto fino a circa 130 cavalli, ma poi conviene prendere una macchina turbo anzichè standard! Questo probabilmente perché con l’aspirazione turbo si risparmia qualcosa in termini di costo dell’ infrastruttura del motore e prestazioni.
Una volta fatta un’analisi esplorativa dei dati, decidiamo di provare a creare un modello di regressione lineare che provi a descrivere in modo significativo i dati che abbiamo a disposizione.
Iniziamo con una regressione lineare semplice, prendendo come variabile covariata Horsepower, che sappiamo avere una correlazione alta col prezzo.
reg<-lm(log(Price)~Horsepower, data=mydata)
summary(reg)
##
## Call:
## lm(formula = log(Price) ~ Horsepower, data = mydata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.64356 -0.18862 -0.06348 0.19537 0.81976
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.1872907 0.0589965 138.78 <2e-16 ***
## Horsepower 0.0112502 0.0005354 21.01 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2816 on 191 degrees of freedom
## Multiple R-squared: 0.698, Adjusted R-squared: 0.6965
## F-statistic: 441.5 on 1 and 191 DF, p-value: < 2.2e-16
Notiamo che le stime dell’intercetta e dello slope della retta di regressione sono significative dal punto statistico, e che il modello è accettabile avendo un R2 elevato.
Dopo una serie di test abbiamo deciso di migliorare il modello aggiungendo altre variabili correlate al prezzo.
linear <- lm(log(Price) ~ Make + Aspiration + BodyStyle + WheelBase + Length + Width +
Height + CurbWeight + NumOfCylinders + EngineSize + PeakRpm +
FuelType + EngineLocation, data = mydata)
summary(linear)
##
## Call:
## lm(formula = log(Price) ~ Make + Aspiration + BodyStyle + WheelBase +
## Length + Width + Height + CurbWeight + NumOfCylinders + EngineSize +
## PeakRpm + FuelType + EngineLocation, data = mydata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.269327 -0.054891 0.002972 0.061271 0.281341
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.696e+00 8.493e-01 9.061 5.74e-16 ***
## Makeaudi 6.723e-02 1.191e-01 0.565 0.573214
## Makebmw 2.983e-01 8.938e-02 3.338 0.001061 **
## Makechevrolet -3.259e-01 1.141e-01 -2.857 0.004877 **
## Makedodge -3.790e-01 9.153e-02 -4.141 5.70e-05 ***
## Makehonda -2.156e-01 9.003e-02 -2.395 0.017832 *
## Makeisuzu -3.933e-01 1.111e-01 -3.540 0.000530 ***
## Makejaguar -4.413e-01 1.303e-01 -3.386 0.000902 ***
## Makemazda -1.350e-01 8.658e-02 -1.559 0.121028
## Makemercedes-benz -1.511e-01 1.425e-01 -1.060 0.290705
## Makemercury -2.059e-01 1.445e-01 -1.425 0.156332
## Makemitsubishi -4.379e-01 8.630e-02 -5.074 1.11e-06 ***
## Makenissan -2.003e-01 8.112e-02 -2.469 0.014660 *
## Makepeugot -4.509e-01 1.059e-01 -4.257 3.61e-05 ***
## Makeplymouth -3.982e-01 9.064e-02 -4.393 2.08e-05 ***
## Makeporsche 2.196e-01 1.468e-01 1.496 0.136638
## Makesaab 3.769e-02 9.866e-02 0.382 0.702940
## Makesubaru -2.967e-01 8.582e-02 -3.457 0.000706 ***
## Maketoyota -2.651e-01 7.940e-02 -3.338 0.001059 **
## Makevolkswagen -1.375e-01 8.710e-02 -1.579 0.116374
## Makevolvo -1.794e-01 9.622e-02 -1.865 0.064161 .
## Aspirationturbo 1.267e-01 3.053e-02 4.149 5.53e-05 ***
## BodyStylehardtop -1.447e-01 6.943e-02 -2.085 0.038748 *
## BodyStylehatchback -1.839e-01 6.337e-02 -2.903 0.004248 **
## BodyStylesedan -1.195e-01 6.653e-02 -1.797 0.074341 .
## BodyStylewagon -1.572e-01 7.470e-02 -2.105 0.036945 *
## WheelBase 2.007e-02 5.229e-03 3.838 0.000181 ***
## Length -4.986e-03 2.955e-03 -1.687 0.093587 .
## Width 1.049e-02 1.329e-02 0.789 0.431355
## Height -3.548e-02 7.607e-03 -4.663 6.73e-06 ***
## CurbWeight 6.897e-04 8.051e-05 8.566 1.08e-14 ***
## NumOfCylindersfive -1.176e-01 1.080e-01 -1.088 0.278100
## NumOfCylindersfour -5.768e-02 1.429e-01 -0.404 0.686998
## NumOfCylinderssix -1.035e-01 1.370e-01 -0.756 0.450855
## NumOfCylindersthree 1.845e-01 2.044e-01 0.903 0.368004
## NumOfCylinderstwelve -3.097e-02 1.930e-01 -0.160 0.872712
## EngineSize 7.446e-04 9.128e-04 0.816 0.415959
## PeakRpm 4.632e-05 2.846e-05 1.627 0.105751
## FuelTypegas 5.321e-02 4.092e-02 1.300 0.195438
## EngineLocationrear 5.460e-01 1.566e-01 3.486 0.000639 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1095 on 153 degrees of freedom
## Multiple R-squared: 0.9634, Adjusted R-squared: 0.9541
## F-statistic: 103.3 on 39 and 153 DF, p-value: < 2.2e-16
Il modello ottenuto è molto soddisfacente è possiamo ritenere che sia significativamente predittivo sulla variabile risposta price.
Per esserne sicuri comunque andiamo a vedere l’analisi dei residui.
plot(linear,pch = 19)
Cerchiamo di capire se i residui soddisfano le proprietà di omschedasticità:
1) Il grafico tra i fitted values (il prezzo) e i residui ha una retta di regressione che è quasi orizzontale. Idealmente vogliamo che la distribuzione dei residui non cambi rispetto ai fitted value (B0=0)
2) Sembra inoltre che i residui si distribuiscano come una normale giudicando il grafico quantile quantile, cioè stanno piu o meno sulla bisettrice.
3) grafico fitted values e residui standardizzati, anche qui la retta di regressione è orizzontale e vediamo che praticamente tutti i dati dati rientrano nell’intervallo \[-1.5,+1,5\] della gaussiana a parte qualche piccolo outlier.
4) Nell’ultimo grafico vediamo che gli outliers stanno dentro alla distanza di Cook. Quindi la loro leva non crea grossi problemi al modello, e possiamo ritenerci soddisfatti.
Vogliamo individuare dei sottogruppi di osservazioni che siano omogenee secondo un determinato criterio.
In particolare vogliamo vedere come possiamo suddividere il dataset in base al prezzo dell’auto e consumo in città, che potrebbe essere una domanda interessante che si potrebbe porre un compratore.
Visualizziamo innanzitutto la distribuzione dei dati e poi cerchiamo una serie di gruppi simili tra loro.
ggplot(data = mydata, aes(x=CityMpg,y=Make))+ geom_point(colour="darkblue", fill="lightblue") + ggtitle("CityMpg vs Make")
ggplot(data = mydata, aes(x=Price,y=Make))+ geom_point(colour="darkblue", fill="lightblue") + ggtitle("CityMpg vs Price")
Decidiamo di utilizzare il metodo delle K-medie, cercando di trovare un’indicazione del K migliore da utilizzare come riferimento per il metodo dei medoidi.
df0<- mydata[, c("Make","CityMpg","Price")]
df1<- mydata[, c("CityMpg","Price")]
crit<-0
for (i in 2:10 ) {
set.seed(7)
mydatagroup<- kmeans(scale(df1),i, nstart = 10)
crit[i-1]<-mydatagroup$tot.withinss
}
plot(1:9,crit, pch=19, xlab = "K", ylab = "whitinss", col="darkblue")
Abbiamo scelto di utilizzare K=4 per il metodo dei medoidi.
pam.out<- pam(df1,4, metric = "euclidean", stand = T)
pam.out$medoids
## CityMpg Price
## 116 19 16630
## 167 26 9538
## 73 16 35056
## 93 31 6849
plot(df1, col=(pam.out$cluster+1), main="PAM {K=4}",pch=20)
points(pam.out$medoids, pch=as.character(pam.out$cluster[pam.out$id.med]))
clusters<-as.factor(pam.out$cluster)
ggplot(data = df0, aes(x=CityMpg,y=Make))+ geom_point(aes(colour=clusters))
ggplot(data = df0, aes(x=Price,y=Make))+ geom_point(aes(colour=clusters))
plot(pam.out, which=2, main="", col = "green")
L’avarage silhouette è soddisfacente per aver scelto K=4.
Abbiamo trovato relazioni interessanti tra le variabili, e abbiamo testato la significatività di queste relazioni.
Abbiamo notato che LogPrice può essere interpretato come una normale.
Abbiamo scovato informazioni interessanti tra le variabili correlate al prezzo.
Abbiamo creato un modello di regressione lineare che possa predire, in base ai dati, l’andamento del prezzo.
Abbiamo creato un suddivisione del dataset in base al consumo di carburante in città, aiutando un possibile compratore nella scelta dell’acquisto dell’auto.